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ABSTRACT 

The  phenomenon  of  thermally  sustained  pressure  oscil- 
lations which  frequently  occur  over  liquid  helium  and 
other  liquified  gases  is  modeled  and  analyzed.   The 
model  is  a  distributed  system  with  a  linear  temperature 
gradient.   The  interaction  between  the  vapor  motion  and 
heat  transfer  is  analyzed  to  derive  an  expression  for  the 
time  history  of  the  pressure  oscillations. 

Cyclic  interpretation  of  the  pressure  history  results 
in  a  relationship  between  four  parameters  which  control 
the  behavior  of  the  oscillations.   The  four  parameters 
contain  information  on  amplitude  of  motion,  slenderness 
ratio  of  the  tube,  characteristic  lengths  of  the  tube, 
thermal  conductivity,  specific  heat  ratio  and  viscosity  of 
the  gas  undergoing  the  oscillations  and  boundary  layer  decay. 

The  effect  of  changes  within  the  parameters  on  the 
theoretical  behavior  of  the  oscillating  system  shows  good 
agreement  with  the  behavior  of  oscillations  in  previous 
experimental  apparatus  when  similiar  changes  are  made. 
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Nomenclature 

The  following  nomenclature  will  be  used  throughout 

the  text  unless  otherwise  noted. 

Symbol  Definition 

A  Cross  sectional  area  of  tube, 

A*  Surface  area  for  heat  transfer. 

Aw  Surface  area  for  shear  stress. 

A-j  Constant  coefficient. 

A2  Constant  Coefficient. 

a  Temperature  gradient  along  tube  wall. 

C  Constant. 

Cp  Specific  heat  at  constant  pressure. 

cv  Specific  heat  at  constant  volume. 

D  Viscous  drag  coefficient. 

D^  Constant  coefficient. 

E>2  Constant  coefficient. 

d  Diameter  of  tube, 

k  Thermal  conductivity. 

L  Length  of  tube. 

L1  Reduced  length  for  shear  stress. 

1  Length  of  temperature  gradient. 

lf  Length  of  compressible  volume, 

m  Mass. 

Nj  Non-dimensional  parameter. 

Kg  Non-dimensional  parameter. 
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Symbol  Definition 

N3  Non-dimensional  parameter. 

n  Angular  frequency. 

P  Pressure. 

P  Average  pressure. 

P«a  Storage  vessel  pressure. 

Q  Heat  transfer  rate. 

R  Gas  constant. 

T  Temperature. 

T]_  Cold  vapor  temperature. 

Th.  Warm  volume  temperature. 

Tw  Wall  surface  temperature. 

T  Average  temperature. 

t  Time. 

U  Internal  energy. 

u  Velocity. 

u0  Maximum  velocity  amplitude. 

V  Volume . 

v  Average  velocity, 

x  Direction  along  tube  or  flat  plate. 

xm  Thullen's  maximum  position  of  inertial  element. 
x0     .  Maximum  displacement  amplitude. 

y  Height  above  flat  plate. 

ym  Maximum  height  above  flat  plate. 

ot  Thermal  diffusivity. 

3*  Specific  heat  ratio. 

Spacial  decay  rate  of  momentum  boundary  layer. 
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Symbol  Definition 

v  Kinematic  viscosity. 

>**•  Viscosity. 

¥  Spacial  decay  rate  of  thermal  boundary  layer. 

J>  Density. 

Xw  Shear  stress  on  tube  wall. 

(p  Non-dimensional  parameter. 

\Jf  Non-dimensional  parameter. 


I.   INTRODUCTION 

Thermally  sustained  pressure  oscillations  are 
often  found  in  the  vapor  trapped  within  cryogenic  equip- 
ment.  Particularly  susceptable  to  this  phenomenon  are 
the  tubes  and  piping  which  penetrate  storage  containers 
of  liquified  gases.   These  tubes  are  characterized  by 
being  closed  at  their  outside  end  which  is  usually  at 
room  temperature  and  open  at  their  other  end  which  is  in 
communication  with  the  cold  vapor  of  the  liquified  gas. 
Since  it  is  usually  necessary  to  fill,  vent,  or  measure 
temperature  and  pressure  within  storage  vessels  by  means 
of  a  tube  or  pipe  penetration,  construction  of  cryogenic 
equipment  without  this  configuration  is  nearly  impossible. 

Occasionally  the  oscillations  have  been  used  to  good 
advantage  for  such  purposes  as  vapor  mixing  or  liquid 
level  measurement.   However,  in  most  ins+ances,  the 
presence  of  the  pressure  oscillations  is  undesirable. 
When  their  pressure  amplitudes  are  sufficiently  large, 
the  oscillations  can  make  vapor  pressure  measurements 
erroneous.   Since  the  open  end  of  the  tube  is  in  free 
communication  with  the  cold  vapor  above  the  liquid,  there 
is  a  heat  pumping  action  to  the  low  temperature  region 
which  greatly  increases  the  boil-off  rate  of  the  liquified 
gas  within  the  container. 

Encounters  with  the  phenomenon  in  cryogenic  apparatus 
is  well  documented  in  cryogenic  engineering  literature. 
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The  phenomenon  is  moot  readily  observable  over  liquid 
helium.   Investigations  by  Ditmars  and  Furukawa  (3)  and 
Thullen  (7)  have  shown  that  the  pressure  oscillations 
may  occur  over  liquid  hydrogen  or  liquid  nitrogen,  although 
their  detection  is  somewhat  more  difficult  than  over 
liquid  helium.   Thullen  (7)  gives  a  rather  complete  survey 
of  encounters  up  to  the  present  time. 

Essentially  since  the  late  1940* s,  the  literature 
has  dealt  with  the  description,  detection,  and  damping  of 
the  oscillations  due  to  their  unv/anted  nature.   Curosiously 
enough,  it  is  Lord  Rayleigh  (5)  that  is  most  often  quoted 
in  explaining  the  oscillations.  In  his  work  The  Theory 
of  Sound ,  Volume  II,  under  the  section  entitled  "Maintenance 
by  Keat  of  Aerial  Vibrations,"  Lord  Rayleigh  describes 
a  form  of  thermally  sustained  vibrations  encountered  by 
glass  blowers  when  blowing  a  bulb  at  the  end  of  a  long 
narrow  tube.   Just  as  found  in  cryogenic  applications, 
the  tube  was  communicating  between  temperature  extremes 
and  the  tube,  itself,  certainly  must  have  possessed  a 
large  temperature  gradient  along  its  wall.   Qualitatively 
Lord  Rayleigh  states  that,  as  the  gas  is  moved  from  the 
cooler  regions  of  the  tube  to  warmer  regions,  "the  adjustment 
of  temperature  takes  time  (italics),  and  thus  the  temp- 
erature  of  the  air  deviates  from  that  of  the  neighboring 
parts  of  the  tube,  inclining  towards  the  temperature  of 
that  part  of  the  tube  from  (italics)  which  the  air  has  just 
come.   From  this  it  follows  that  at  the  phase  of  greatest 


condensation  heat  is  received  by  the  air  and  at  the 
phase  of  greatest  rarefaction  heat  is  given  up  from  it, 
and  thus  there  is  a  tendency  to  maintain  the  vibrations." 
Lord  Rayleigh  further  points  out  that  the  bulb  and  the 
great  range  of  temperature  is  necessary  for  the  maintenance 
of  the  vibration.   Within  cryogenic  apparatus  construction, 
low  conductivity  material  such  as  glass,  cupro-nickel, 
or  stainless  steels  are  used  for  the  tubes.   Thus,  a 
severe  temperature  gradient  must  exist  near  the  cold  end 
of  the  tube;  and  the  warm  end  of  the  tube  and  associated 
valving  act  as  the  bulb  of  Lord  Rayleigh1 s  high  temperature 
phenomenon. 

Extensive  experimentation  has  been  conducted  to 
determine  the  types  of  materials  and  geometric  configurat- 
ions which  will  foster  and  sustain  the  pressure  oscillat- 
ions in  a  room  temperature  to  cryogenic  temperature  environ- 
ment.  Clement  and  Gaffney  (2),  Ditmars  and  Furukawa  (3), 
Bannister  (1),  and  Thullen  (7)  have  in  their  own  ways 
measured  frequencies,  pressure  amplitudes,  and  devised 
methods  of  damping  the  oscillations.  From  their  experiment- 
ions  and  reports  it  is  evident  that  the  tube  diameter, 
length,  and  magnitude  of  the  longitudinal  temperature 
gradient  in  the  tube  wall  are  the  most  important  physical 
characteristics  of  the  apparatus. 

Analytical  investigations  have  yet  to  lead  to  a 
satisfactory  theoretical  explaination  of  the  phenomenon 
associated  with  cryogenic  apparatus.   Norton  and  Kulenhaupt 
(4)  have  produced  a  computer  program  to  simulate  the 


pressure  oscillations.   In  their  program  the  continuous 
temperature  gradient  was  approximated  by  a  large  number  of 
constant  temperature  regions.   The  program  gave  encouraging 
results  when  compared  to  experimental  measurements  of  the 
pressure  oscillations.   Thullen  (7)  proposed  a  lumped 
parameter  model  in  which  he  was  able  to  form  two  dimension- 
less  variables,  one  containing  frequency,  friction  and 
temperature  gradient  information,  and  the  other  containing 
amplitude  and  heat  transfer  information. .  From  this 
approach  the  exact  behavior  of  any  distributed  system 
can  not  be  predicted;  but,  it  does  point  out  the  controlling 
variables.   Both  the  approaches  made  by  Norton  and  Mulen- 
haupt  and  Thullen  have  added  significantly  to  the  theo- 
retical understanding  of  thermally  sustained  pressure 
oscillations  in  cryogenic  equipment.   However,  due  to  the 
approaches  taken,  neither  theory  was  able  to  capture  the 
inherent  distributiveness  of  the  system  and  the  continuous 
nature  of  the  temperature  gradient  along  the  tube  wall. 

It  Is  the  intention  of  this  author  to  define  a  model 
for  the  pressure  oscillations  which  encompasses  the  im- 
portant points  of  previous  experimentation,  and  analyze 
the  model  *n  such  a  manner  as  to  retain  as  much  distribut- 
iveness of  the  system  and  continuity  of  the  temperature 
gradient  as  possible,  while  maintaining  a  manageable 
theoretical  discription  of  the  phenomenon. 


II.   ANALYSIS 

1.   The  Model. 

The  model  chosen  for  the  theoretical  analysis  of 
the  thermally  sustained  pressure  oscillations  closely 
approximates  the  geometric  configuration  often  described 
in  the  literature.   Common  points  in  these  configurations 
are  the  closed  volume  at  the  outside  end  of  the  pene- 
trating tube,  the  longitudinal  temperature  gradient  along 
the  tube  wall,  and  the  free  communication  of  the  cold 
vapor  with  the  open  end  of  the  tube.   The  system  for 
analysis  is  illustrated  in  Figure  1. 
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Figure    1.      Model  For  Analysis 


It  is  Important  that  the  temperature  gradient 
along  the  tube  wall  be  of  sufficient  strength  to  allow 
heating  of  the  vapor  upon  compression  and  cooling  upon 
expansion.   The  maximum  amplitude  of  the  oscillations 
is  limited  due  to  viscous  friction.   Thus  the  model  of 
real  systems  has  all  the  elements  to  foster  the  growth  of, 
and  sustain,  pressure  oscillations.   Also,  the  distributed 
nature  of  the  system  and  the  continuity  of  the  temperature 
gradient  have  been  preserved. 

Within  cryogenic  apparatus  the  oscillations  are 
spontaneous  in  origin;  and,  after  a  period  of  time,  their 
amplitudes  grow  to  a  maximum  value.   The  analysis  of  this 
model  will  assume  that  a  sinusoidally  varing  steady  state 
condition  has  been  reached  within  the  system.   The  re- 
lationships between  time,  motion,  and  heat  transfer  will 
be  combined  to  define  the  pressure  as  a  function  of  time, 
physical  characteristics  of  the  vapor,  and  geometric 
configuration. 

2»   Derivation  of  the  Time  and  Pressure  Relationship. 

When  the  model  is  considered  as  an  open  system, 
the  oscillations  of  the  vapor  within  the  tube  cause 
energy  to  be  added  or  lost  to  the  system  by  virtue  of 
mass  entering  or  exiting  the  open  end  of  the  tube,  and 
by  heat  transfered  to  the  vapor  because  of  the  temperature 
gradient  along  the  tube  wall.   The  first  law  of  thermo- 
dynamics for  an  open  system  applies  to  the  control  volume 


within  the  closed  tube: 


dt   cPTldt  +  Q  (1) 


The  time  rate  of  change  of  mass  entering  the  system 
across  the  open  end  of  the  tube  is  given  by: 


a!=/Av(t>  (2) 


where  v(t)  is  the  average  velocity  of  the  fluid  particles 
at  the  end  of  the  tube  at  a  given  time  t. 

The  perfect  gas  relations  and  the  definition  of 
internal  energy  give  the  rate  of  change  of  pressure  in 
terms  of  the  rate  of  change  of  internal  energy: 

(3) 

W 

(5) 

P  =  CU  (6) 


p 

= 

mRT 
V 

u 

a 

mcvT 

p 

= 

-5-u 

OyV 

where  C  = 


_  R 


cvV 


Since  both  pressure  and  internal  energy  are  functions 
of  time : 


§  -  "8  ™ 
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When  equation  (7)  is  substituted  into  equation  (1), 
the  rate  of  change  of  pressure  is  given  as  a  function  of 
time : 


d? 

at 


=  C  fcp^TivCt)  +  i(t)]  (8) 


Density  is  not  a  variable  of  primary  interest;  and, 
over  a  cycle  of  the  oscillations,  its  variation  is  small. 
Therefore,  the  following  simplifing  assumption  and  sub- 
stitution is  made  in  equation  (8): 

f  =  h  (9) 

RT 

where  P  is  the  average  pressure  at  an  average  temperature, 
T,  during  a  cycle : 


dP  =  cf 
dt 


-^-(^1)v(t)  +  Q(t) 


(10) 


When  v(t)  and  Q,(t)  are  known,  equation  (10)  can  be 
integrated  to  give  the  pressure  as  a  function  of  time. 

Since  the  usual  tube  diameters  are  large  in  comparison 
to  the  thermal  and  momentum  boundary  layer  thicknesses, 
the  derivation  of  the  velocity  term,  v(t),  and  the  heat 
transfer  term,  Q,(t),  in  equation  (10)  will  be  carried  out 
in  a  rectangular  co-ordinate  system.   Appendix  A  shows  the 
range  of  tube  diameters  for  which  this  approach  will  be 
valid. 


3.   Solution  of  the  Velocity  Term,  v(t). 

Schlichting  (6)  presents  the  boundary  layer  solution 
for  the  flow  of  a  fluid  near  an  oscillating  flat  plate. 
For  the  purpose  of  this  analysis,  Schlichting' s  solution 
will  be  modified  to  fit  the  boundary  condition  that  the 
fluid  is  sinusoidally  oscillating  over  a  stationary  flat 
plate.   Accordingly,  the  velocity  distribution  in  the 
fluid  above  the  flat  plate  as  a  function  of  height  above 
the  plate,  y,  and  time,  t,  is  given  by: 

u(y,t)  =  u0e~*ycos  (nt  -  <Ty)  -  u0cos  nt       (11) 

where  £  =V— H  *s  the  spacial  decay  rate  of  the  momentum 
boundary  layer. 

The  average  velocity  of  the  fluid  entering  the  tube 
is  obtained  by  assuming  that  the  average  velocity  of  the 
fluid  over  the  flat  plate  to  a  given  distance  is  analogous 
to  that  entering  the  tube: 

ym 


r  Jm 

t)  =  Jou(y>t^dy 


v(t)  =  jo»^»wuy  (12) 

ym 

where  ym  corresponds  to  the  tube  radius,  d/2,  and  is 
sufficiently  large  so  that  in  equation  (1?)  e"Jym<$^1. 

When  the  average  velocity  is  determined  in  this  manner, 
the  final  solution  of  the  velocity  term  is: 


v(t)  =  u, 


(  — A-  -  1)cos  nt  +    1   sin  nt 

ym 


^m  2<T 


(13) 
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4.   Solution  of  the  Heat  Transfer  Term.  £. 

The  model  assumes  that  the  only  heat  transfer  to  and 
from  the  vapor  Is  In  the  area  of  the  tube  where  the 
temperature  gradient  exists.   The  rest  of  the  volume  of 
the  tube  is  assumed  to  be  adiabatic  with  respect  to  the 
vapor.  As  with  the  velocity  term  derivation,  the  flat 
plate  approach  is  assumed  to  be  applicable. 

To  begin,  an  energy  balance  is  made  on  a  control 
volume  of  unit  depth  positioned  above  the  flat  plate. 
There  is  a  vertical  heat  flux  through  the  control  volume 
and  a  horizontal  enthalphy  flux.   Figure  2  illustrates 
the  energy  balance. 


j>  udyCpT 


y 


-kdx|l  +        *   (~kdx|l)dy 

©y      A^y       «y 


r 


Unit  Depth 
Control  Volume 


-s* 


L 


-> 


j3udycDT  + 

^r(/>udycpT)dx 


I 


t  -kdxil 

77777777777777777777777 


-e-x 


Figure  2.   Energy  Balance  For  Heat  Transfer  Derivation 
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The  energy  balance  on  the  control  volume  yields  the 
following  differential  equation: 

^<4i  "  4?*  w<y.*l  (14) 

Equation  (1A)  may  be  rearranged  into  a  diffusion  equation 
with  a  velocity  type  forcing  function: 

I  $1   =  e*Jx  -  u(y,t)^T  (15) 

2^t      ay2         *x 

c_         . 
where  3*  -   -i-  and  c^<  =  £— . 

cv         /cp 

When  equation  (15)  has  been  solved  for  the  temperature 
distribution  in  the  fluid  above  the  flat  plate,  the 
temperature  gradient  at  the  wall  can  be  found  by  differ- 
entiating the  solution  with  respect  to  y  and  setting  y 
equal  to  zero.   This  will  enable  the  heat  transfer  term 
to  be  realized  as  a  function  of  time. 

The  following  assumptions  are  made  in  reguard  to  the 
solution  of  equation  (16): 

a,  A  linear  temperature  gradient  exists  in  the  wall 
so  that  the  wall  surface  temperature  is  given  by: 

Tw  =  ax  +  E  (16) 

b.  The  same  linear  temperature  gradient  as  in  the 
wall  exists  within  the  gas  so  that  the  form 

of  the  temperature  distribution  in  the  gas  above 
a  point  on  the  plate  is  given  by: 

T(y,t)  =  f(y,t)  +  ax  (17) 
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c.  There  is  neglible  heat  conduction  in  the  gas 
in  the  x-direction. 

d.  The  velocity  distribution  is  given  by  equation 
(11)  and  the  average  velocity  is  given  by 
equation  (13). 

Equation  (17)  must  also  fit  the  boundary  condition 
that  at  the  wall  the  wall  surface  temperature  and  the 
gas  temperature  must  be  equal. 

When  equation  (17)  is  substituted  into  equation  (15), 
the  differential  equation  which  determines  the  temperature 
distribution  above  the  plate  is: 

£sC=e4§2-au(Jr't)  (18) 

where  the  function  f(y,t)  is  subject  to  the  boundary 
condition  that  at  y  =  0,  f(0,t)  =  E. 

The  solution  equation  (18)  is  of  the  form: 

f(y,t)  =  e~*y  AjSln  (nt  -  <*y)  +  A2sin  (nt  -  Sj) 

+  B,sin  nt  +  13  cos  nt 
1  2 

+  e~*y  D^infnt  -  Ty)  +  D2cos(nt  -/y) 
+  constant  (19) 

Appendix  B  shows  the  determination  of  constants  associated 
with  equation  (19). 

The  complete  solution  to  equation  (17)  giving  the 
temperature  distribution  in  the  gas  above  the  flat  plate  is: 
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T(y,t)   =  E   +  ax   +  u0a  e"   ^     sin(nt   -  &y) 

2*£i  -  n 

+  u0aj£sin  nt 
n 

-  2u0ae<($2^7n  e~^sin(nt  -  fj)  (20) 

2*tJz-   n 
? 

where  f  =-\[    n   ,  soacial  decay  rate  of  thermal  boundary  layer 

2**r . 

S  =V  n"~  ,  spacial  decay  rate  of  momentum  boundary 
2  "V  layer. 

The  heat  transfered  from  the  plate  to  the  oscillating 
vapor  is  given  by: 

Q  =  -k^Tl  A'   *  (21) 

fcyjy  =o 

When  equation  (20)  is  differentiated  with  respect  to 
y  and  evaluated  at  y  =  0  and  the  result  substituted;  into 
equation  (21 ),  the  heat  transfer  term,  Q,  is  given  by: 

Q  =  -kaA'uptf  2^Sf^  -    1  [sin  nt  +  cos  nt       (22) 


5.   Final  Solution  of  the  Pressure-Time  Relationship. 
With  the  results  of  sections  3.  and  4.  combined 
with  the  differential  equation  of  section  2.,  the  final 
derivation  of  the  gas  pressure  as  a  function  of  time  can 
be  determined. 
Equation  (10)  of  section  2.  is: 
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dP 

dt 


-   C 


cLpPA 


T 


v(t)    +  Q,(t) 


Equation    (13)    of  section  3.    is: 


v(t)   =  u0 


J ilcos   nt   +  _1    sin  nt 


Equation    (22)    of  section  4.    is: 


2*4  «r 


A&o     h*^?-    il 
r«  -  ».  L  n         J  L 


sin   nt  +  c 


os  nt 


When  equations  (13)  and  (22)  are  substituted  into  (10) 
and  equation  (10)  is  integrated  with  respect  to  time,  the 
gas  pressure  as  a  function  of  time  is: 


P(t)  =  C 


R 
kaA 


MM 


+  P. 


sin  nt  -   1   cos  nt 

2*yffi 

1  l-cos  nt  +  sin  ntl 


(23) 


where  Poo  is  the  constant  of  integration  or  pressure  at 
time  t  =  0". 

With  the  following  definitions  the  compressible 
volume  prersure  given  by  equation  (24)  can  be  made  non- 
dimensional  and  more  manageable  for  practical  applications: 
Let      N«  = 


0   2^7m 
N1  =  2ctSfg 


n 


N2  =  2*&£ 

n 
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t  =  a. 


H  =  2*^ 


Upon  the  substitution  of  the  above  definitions, 
equation  (23)  becomes: 


P(t) 
UBCUq/iS 


■[ 


HI* 
B 

-  |ht*. 

LB 


(N0  - 


No 


'    N2  -  1 


sin  nt 


K4  -  1 


N2  -  1 


cos  nt 


(24) 


In  a  likewise  manner  the  average  velocity,  v(t),  can 
be  made  non-dimensional: 


Z&1  ~    (N0  -    i)cos  nt  +  N0sin  nt 

UQ 


(25) 


III.   DISCUSSION 

With  the  pressure  history  of  the  warm  volume  and 
the  motion  of  the  mass  taking  part  in  the  oscillations 
now  known,  a  dynamical  approach  to  the  system  will  be 
used  to  derive  the  relationship  between  the  parameters 
which  control  the  oscillations.   The  model  of  the  system 
for  this  investigation  will  be  similar  to  that  used  by 
Thullen  (7)  and  is  shown  in  Figure  3. 


— >j  d  k — 


Viscous 

Damping 


it 


x.v 


Pp(t)    Temperature 
Profile 


■**2TW 


Shear  Stress 
on  Tube  Wall 


Figure  3.   Model  For  Cyclic  Interpretation 
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Newton's  law  of  motion  applies  to  the  inertial 
element  of  the  model : 

mx  +  Dx  -  A[p(t)  -  *e(tj\    =   0  (26) 

The  transient  solution  to  equation  (?6)  is  not  of 
Interest  and  can  be  avoided  by  the  cyclic  integration 
of  equation  (26)  with  respect  to  displacement.   The 
net  work  done  in  accelerating  the  inertial  element  through 
a  cycle  is  zero: 


/ 


mxdx  =  0  (27) 


The  two  remaining  terms  of  equation  (26)  are  now 
equated  and  point  out  that  the  cyclic  work  done  on  the 
volume  to  the  left  of  the  inertial  mass  is  equal  to  the 
cyclic  work  done  against  friction: 

OA[p(t)  -  Pe(t)]dx  =  6Dxdx  (28) 

The  identity  dx  =  <*2£  dt  =  v(t)dt  is  substituted 

dt 
into  equation  (28)  so  that  the  cyclic  integration  can 

be  performed  with  respect  to  time  instead  of  displacement. 

i  A[p(t)  -  Pe(t)]v(t)dt  =^Dv2(t)dt  (29) 

If  the  exit  pressure  is  constant,  Pe(t)  =  P<*>,  the 
following  enequallty  must  be  satisfied  for  the  oscillations 
to  be  self  sustaining: 

OAfp(t)  -  pJ  v(t)dt=0Dv2(t)dt  (30) 
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In  substituting  equations  (24)  and  (25)  into  (30) 
and  performing  the  cyclic  integration,  the  following 
relationship  is  established: 


ABC 


Nj  -  1 


*(N0  -  1)2  +  BJ 


nD  No  -  1 

Thullen's  (7)  non-dimensional  parameter  containing 
frequency,  friction,  and  temperature  gradient  information 
is  as  follows: 


(31) 


N.  =  ?^2#   b  _i_ 


8      Dt^TT 


(32) 


Upon  substitution  of  the  following  definitions 
of  the  variables  in  equatiion  (32),  the  Thullen  parameter 
becomes : 

P0  =  /RTi 
A  =  TTd^ 


s  =  ira?V 
"T" 


b  =  AT 
T. 

u>  =  n 


N, 


d< 


P 


RaT  a* 


Dn 


(33) 
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Upon  sutsitituion  of  the  actual  variables  for  B  and 

C  in  the  non-dimensional  scrou^ing  ABC  in  equation  (31 )  and 

nD 
makins  the  following  definitions,  this  non-dimensional 

group  is  seen  to  contain  Thullen's  parameter  plus  additional 

information: 


N* 

=  ABC 
nD 

— 

AkaA'cf  3" 

n'^D 

R 
cvV 

A 

=  7fd2 

4 

A' 

=  TTdl 

V 

=  7Td2 

~5~ 

V 

i 

cT   =V    n  =~\lnf "=  P"\Tn 
20       2^        y*p 


(34) 


N*  = 


d2 


j>r^t^ 


Dn 


4iTk 
c^nf  &iJ2/Lf 


(35) 


N*  =  NT 


4-rrk   __  < 
cvnid  -^f2yu/ 


(36) 


In  order  to  carryout  the  cyclic  integration  of 
equation  (2?)  with  a  variable  exit  pressure,  Pe(t), 
further  definition  .of  the  exit  pressure  term  and  the 
damping  coefficient  term,  D,  are  required. 

First,  the  exit  pressure  is  examined: 


C 


?(t) — =&\53i&—  p.(t) 

e-     +  V(t) 

Outflow:  Pe(t)  =  Po 

Inflow:  Pe(t)  =  Pso  -  i/v2(t) 
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<£  APe(t)v(t)dt  =    |APe(t)v(t)dt 
J  -/inflow  |  cycle 

+  JAPe(t)v(t)at  (37) 

*^out   flow  i  cycle 

The  outflow  half-cycle  contributes  nothing  to  the 
cyclic  integration.   In  the  inflow  half-cycle  the  sign 
of  the  velocity  is  negative.   Therefore,  the  cyclic 
integration  of  equation  (29)  reduces  to: 

A  A[P(t)  -  pjv(t)dt  -  r|-Aj>v3(t)dt  =<£r>v2(t)dt     (38) 
J  -'inflow  J  cycle' 

Second,  the  damping  coefficient  term  can  be  determined 
through  a  pipe  flow  analogy  involving  shear  stress: 


r*ii      1^ 


F  =  (P,  -  P2)  A  =TwAw 

^dyjy  =  o 

The  rate  of  change  of  velocity  with  respect  to  y 
is  found  from  the  velocity  equation  (11) 

~l    „  =  -Juofsin  nt  +cob  ntl  (39) 

dy|  y  =  0      °l  J 

The  damping  coefficient  is  now  defined  as: 


A v<iuu q\ sin  nt  +  cos  ntj 


D  =  -F  =  Awcj#u0[sin  nt  +  cos  ntj  (4o) 

^TtT  ~vTtT~ 


The  right  half  of  equation  (38)  now  requires  the 
cyclic  integration  of  the  following  term: 


/ 


Ajfeuoisln  nt  +  cos_nt| v2(t)dt  (*l) 

v(t) 
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With  the  above  substitutions  for  the  exit  pressure 
and  damping  coefficient  terms  and  subsequent  cyclic 
integration  of  equation  (33),  the  relationship  bet\^een 
the  parameters  which  contribute  to  the  sustaining  of  the 
pressure  oscillations  can  be  established: 

ABC  ug-TrfNi  -  l]  _  2  .Mug  r,N    n2   N2  ]  3/2  _ 

Hft«JLa§  [(M-0  -  1)2  +  N^]1  2      (42) 

Upon  substitution  of  the  following  definitions 
and  previous  values  of  B  and  C,  equation  (42)  takes 
the  form: 


*  ■fe-^l- 


fl/??*-  1 

T&-   1 


u0  =  nx0 

A     =  tTdL1  ,    where  L*    is   the   reduced  length    (L  -  i-11) 

W 

of  the  tube  since  the  shear  stress  Is 
not  uniform  over  the  entire  length  of 
the  tube  of  Figure  3. 

2-*   V  2/^ 

H^-MSfe-'W**  -f"    (43) 
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Figure  4.  shows  the  non-dimensional  displacement 

amplitude,  1   2.V—  xft,  slotted  against  the  driving  force 
3  I!  2^ 

oarameter,  j^ATy(g*  -  1] with  lines  of  constant   (//"  . 

l'L>n2  T 

The  graphical  representation  of  equation  (43)  gives 

indication  of  a  given  system's  performance  and  how  that 

performance  will  change  if  the  various  parameters  are 

altered.   The  manner  of  change  suggested  follows  closely 

that  which  has  so  often  been  described  in  literature  on 

thermally  sustained  pressure  oscillations.   Due  to  the 

lower  limit  on  tube  diameter  imposed  by  the  theoretical 

approach  taken,  the  lines  of  constant  if/   in  Figure  4. 

are  in  aggrement  with  the  theory  for  y    =  1 .   Values  of 

\lf  %   1  are  shown  in  dashed  lines  for  completeness. 

Of  particular  interest  is  the  linear  relationship 

between  the  displacement  amplitude  and  the  slenderness 

ratio,  L'/d.   Bannister  (1)  noted  this  in  his  experimentation 

with  oscillations  over  liquid  helium.   He  concluded  that 

the  pressure  amplitudes  are  directly  proportional  to  the 

slenderness  ratio  of  the  gas  column.   The  Bannister  data 

was  used  to  calculate  values  of  iff  and  the  driving  force 

parameter.   The  calculated  operating  points  are  plotted 

on  Figure  4.   Appendix  C  tabulates  the  constants  and  data 

used.   Figure  5«   shows  the  calculated  amplitude  parameter 

plotted  against  the  measured  maximum  pressure  amplitude. 

The  actual  slenderness  ratio  of  Bannister's  tubes  is  also 

plotted  against  the  maximum  pressure  amplitude  to  show 
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the  linearity  comparison. 
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Figure  5.   Slenderness  Ratio  Comparison 
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Thullen  (7)  also  found  a  non-dimensional  displacement 
amplitude  parameter  in  his  analysis:  xm/21.   The  char- 
acteristic length,  1,  was  the  distance  between  two  ideal 
heat  exchangers  in  his  model;  and,  it  was  viewed  as  the 
distance  traveled  by  the  gas  after  a  flow  reversal  before 
significant  heat  transfer  begins.   However,  the  length 
was  not  well  defined  for  a  distributed  system.   The 
results  of  the  present  analysis  indicate  that  Thullen' s 
"1"  might  have  the  following;  value  for  a  distributed  system: 


i  =  3  vm 

'   d   n 
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IV.   CONCLUSIONS 

A  greater  understanding  of  thermally  sustained 
pressure  oscillations  in  liquified  gas  apparatus  is 
possible  through  a  study  of  vapor  motion  and  heat  transfer 
interaction.   From  such  a  study  a  time  dependent  expression 
for  the  pressure  history  of  a  vapor  in  a  tube  can  be 
derived.   With  cyclic  interpretation  the  relationship 
among  the  control ing  parameters  of  the  system  can  be 
established. 

When  changes  are  made  in  the  parameters,  the  behavior 
of  the  theoretical  relationships  appear  to  follow  closely 
that  which  occurs  when  simlliar  changes  are  made  in 
experimental  apparatus. 

The  rectangular  co-ordinate  system  from  which  the 
motion  and  heat  transfer  interaction  was  derived  places 
a  lower  bound  on  the  diameter  of  the  tube  for  which  the 
theory  applies.   Therefore,  an  analysis  with  the  same 
general  approach  to  the  motion  and  heat  transfer,  but  in 
cylinderical  co-ordinates,  would  be  of  value.   It  is  not 
anticipated  that  this  would  alter  the  basic  relationship 
between  parameters  controling  the  oscillations;  however, 
the  cylinderical  system  might  provide  more  exact  numerical 
constants.   Appendix  D  suggests  an  approach  to  the  analysis 
of  the  phenomenon  of  thermally  sustained  pressure  oscillations 
in  a  cylinderical  co-ordinate  system. 


APPENDIX 
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APPENDIX  A 

In  order  to  show  that  the  flat  plate  analogy  is 
valid  for  the  results  desired,  it  is  sufficient  to  show 
that  the  spacial  decay  rate  is  small  in  comparison  to 
the  tube  radius. 


Let  <fy     =  y/x 

n 

5  x  10"^  ft2/sec 
50  rad/sec 
then   A  =  0.0535  in 


,-yA 


1.00 


1  X  2X  3^4 


Hence,  tubes  with  radii > 0.125  in.  will  be  in  the 
range  for  which  the  theoretical  approach  is  valid. 
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APPENDIX  B 

The  determination  of  the  constants  associated  with 
equation  (19)  which  is  the  solution  to  the  differential 
equation  (18)  is  as  follows: 

Equation  (18):  I  £f    .^ff      ,   ,. 

where  u(y,t)  =  u  e~  ^cos(nt  -  Jy)  -  uQcos  nt 

and  f(y,t)  must  satisfy  the  "boundary  condition:  f(0,t)  =  E. 
Equation  ( 19) : 

f  (y,t)  =  e"^y[AlSin(nt  -  Jy)  +  A2cos(nt  -  (fy)] 

+  B^sin  nt  +  BoCos  nt 

+  e"^DlSin(nt  -  ?y)  +  D2cos(nt  -  f  y)] 

+  constant. 
Thus  the  form  of  the  function  f  is:  f(y,t)  =  f -  +  fg  +  Constant 

Seek  the  solution  of  equation  (18)  by  f . : 

—  1   =  e"   ^lA-jncosCnt   -    S y)    -  A2nsin(nt  -    <fy)l 


-»■  Buncos  nt  +  Bpnsin  nt. 


|-|l   =  2<fV*y  [-A2sin(nt  -     Sy)   +  A^osCnt   -     <?y)l 


Upon  substitution  of  these  partial  derivatives  into 
equation  (18)  and  with  like  coefficients  of  the  sine  and 
cosine   terms   equated: 

A.    =       __upa„_v 
1  (2ct<f2-  n) 
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AP(?c^cr2-  £)  =  0,     k0   =  0 


B,  =  u0a-n- 

B2  =  0 

fi(y,t)  =  Hoae"  ysln(nt  -  <4V)  +  u0a2*cos  nt 

However,  fi, alone,  will  not  fit  the  boundary  condition. 
Hence,  an  additional  function  must  be  added  to  f^  and 
that  function  must  be  a  solution  to  the  homogeneous  part 
of  equation  ( 18)  :  J.  ^f  _  j  Sf 

*  at  -  ay2 

f2(y,t)   =  e"ry[DlSin(nt   -  f y)    +  D2cos(nt  -  f y)] 
+  constant 

~2  -  ne-*V[-D2sin(nt  -  fy)    +  Dtcos(nt   -  fy)] 

|~2  =  2fV  y[-D2sin(nt  -fy)   +  D^osfnt  -fy)] 

Upon  the  substitution  of  the  partial  derivatives 
into  the  homogenous  differential  equation  and  with  like 
coefficients  equated,  f2  becomes: 

D2(2*f2-  £)  =  0      D1(2^fE-|)  =  0 

If  D]  /  0  and  D2  ^  0,  then  (2<*T*-  £)  =  0 
Therefore,      . , 

With  the  boundary  condition  now  applied  to  f(y,t), 
the  coeffients  D^ ,  Dp  and  the  "constant"  of  f2  are  determined: 


30 


constant   + 


**n 


2u0eAd    5 


2*cf2-  R 


Therefore, 

constant  =  E 


7/^ 


+    D. 


i sin   nt   +   D2cos   nt  =  E 


D2  =  0 

f2(y,t)   =  e-ry/-2u  a«fd*-\ 

__n 

2tfcP-  n    " 


sin  nt  +  E 


APPENDIX  C 

The  following  constants  for  helium  gas  (-200°F) 
were  used  in  the  Eannister  data  fit  of  equation  (43): 

c*    =   5.7  x  10~4  ft2/sec 

JUL    =  8.3  x  10~6  Ibm/ft-sec 

^  =  4.0  x  10"4  ft2/sec 

k  =  2.2  x  10~5  ETU/sec-ft-°F 
t   =  1.6? 

Data  from  Bannister's  experimentation: 

n  (cps)   d  (cm)  L  (cm)   P  (atm  x  103)  AT  (°F) 

X   33.0      .287  149      56  480 

O   33.4      .482  149      28  480 

D   33.4      .605  149      18  480 

A   22.7     1.240  213      12  480 

Estimated  from  apparatus  geometric  configuration: 

L'  (cm)   1'  (cm) 

X   105       89 

O  105      89 

a     1 05      89 

A  105      89 
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APPENDIX  D 

The  cylinderical  co-ordinate  analysis  will  proceed 
with  analogous  steps  to  those  performed  in  rectangular 
co-ordinates.   Figure  5.  illustrates  the  system  for  analysis. 
Schlichtlng1 s  solution  for  the  velocity  distribution  of 
an  oscillating  flow  through  a  pipe  is: 

Navier-Stokes  Equation: 

Assume  that  the  pressure  gradient  is  caused  by  a  harmonically 
moving  piston: 

"4  Jr-  =  Kcos  nt 
/»  ox 

Use  complex  notation  and  attribute  physical  significance 

only  to  the  real  part: 

_  1 12  _  Keint 
J>   Sic  ~  *e 

Assume  the  velocity  function  has  the  form: 


u(r,t)  =  f(r)elnt 


R 


Tw  I 


Figure  6.   Cylinderical  Co-ordinate  Model 
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The   solution   to   the   Kavier-Stokes   equation  which 
gives   the   velocity  distribution  within   the   pipe    is: 


u(r,t)  =  -ilielnt[l   -  2ftfeil| 

n      L     JoM)J 


where  i   =V^ 

The  temperature  distribution  in  the  pipe  is  assumed 
to  have  the  following  form: 

T(r,t)  =  f(r,t)  +  ax 

Boundary  condition  for  the  temperarure  distribution  is: 
T(R,t)  =  Tw  =  E  +  ax 

f (R,t)  =  E 

The  differential  equation  resulting  from  energy  considerations 

is: 

1  ^f^T 


V  bt 


r  Brl*  br  J 


-  au(r,t) 


The  ultimate  differential  eouation  to  be  solved  is: 


1  £1 
V  bt 


=  of 


111     St. 

r  br       'br- 


-  au(r,t) 


Assume  the  solution  of  the  differential  equation  if  of  the 
form: 


f(r,t)   =  f,(r,t)    +   f2(r,t)   +  C 

1 


where,  .    . 

f^r.t)   =  -iAelnt 


,  £oMl1 

Jo(Rrf)J 


f2(r,t)   =  -iBe 


int 


Jo(Rf) 


+  C,    where    f  =V-ln 
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The  solution  is: 


^(r.t)  =  - 


f2(r,t)  = 


aK 


1 


aK  <*S*n 


1  - 


1  - 


Jo(RcJ)J 

Jo(rf)] 
Jo(Rf)J 


int 


int 


+  E 


The  temperature  distribution  within  the  pipe  is: 


T(r,t)  = 


-  il   .  J0(rJ)1  +  ee 
A  L    Jo(R<T)J   A 


1  -  Jo(rT) 


J0(Rf)J 


int 


-♦■  E  +  ax  ,  where  jS  =  *K  ,  A  =  -i-^J- H,  &=*£? 

n        ar     n 

The  heat  transfer  from  the  wall  surface  to  the  vapor 
in  the  pipe  is: 


Q  ss.-k  27Trl 


*t1 

5rr  =  R 


Therefore,  the  heat  transfer  term  sought  is: 


Q  =-2tTklR  ~ 
A 


JJiMl  +gf£ijgl)1 

Jo(R«5)    Jo(Rf)J 


in- 


is: 


The  average  velocity  of  the  vapor  entering  the  pipe 
R 


t)  =J° 


v(t) 


v(t)=-i£ 


2tfru(r,t)dr 
ffR2 


2  Jo(R/)1  ^int 


< -  as  $$]  •■ 


Now  both  the  velocity  term  and  heat  transfer  terms 
are  functions  of  time  and  can  be  substituted  into  equation 
(10)  of  section  II. 2, : 


g  =  c[&#ftv(t)+act> 
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Upon  substitution  and  integration  with  respect  to  time, 
and  attaching  significance  to  the  real  part  of  the  result, 
the  pressure  as  a  function  of  time  is: 


P(t)    - 


fcHfl] cos  * 


-  2  klR 


aK    I-  °S&Z)  \   Ui(RT) 


f  J0(Rf)  J 


sin  nt 


The  average  velocity  has  the  following  form  when 
returned  to  real  notation: 


v(t)  = 


£  1*1  -  _£  Ji(R*)l  <   + 

»  L      B  JoTrITJ  sln  nt 


The  above  pressure  history  of  the  compressible  volume 
and  the  average  velocity  has  gouping  of  parameters  which 
appear  similar  to  those  found  in  the  rectangular  co- 
ordinate analysis. 
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